PsL Monthly 1993 December

home *** CD-ROM | disk | FTP | other *** search

/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / pascal / p_mat.exe / PMAT.DOC < prev next >

Wrap

Text File | 1993-01-24 | 48.1 KB | 1,311 lines

P-Mat v1.2 An Turbo Pascal program for Recursive Matrix Algebra Mark Von Tress, Ph.D. PO Box 171173 Arlington TX 76003 Compuserve User ID : 71530,1170 Date: January 30, 1993 You are responsible for what you do with the code. Here is a formal disclaimer: DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING FROM USE OF THIS PROGRAM. Copyright (c) Mark Von Tress 1993 Mat-P 3 Table of Contents I. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . 4 II. FINANCIAL CONDITIONS OF USE . . . . . . . . . . . . . . . 4 III. BASIC ALGORITHMS . . . . . . . . . . . . . . . . . . . . 5 A. Structures and Allocation . . . . . . . . . . . . . . 5 B. The global stack: dispatch . . . . . . . . . . . . . 7 C. Recursion . . . . . . . . . . . . . . . . . . . . . . 7 D. Matrix Assignment . . . . . . . . . . . . . . . . . . 8 E. Parameter Passing . . . . . . . . . . . . . . . . . . 9 IV. FUNCTIONS . . . . . . . . . . . . . . . . . . . . . . . . 10 A. Error Handling . . . . . . . . . . . . . . . . . . . 10 B. Array Allocation or Deallocation . . . . . . . . . . 11 C. Stack Control . . . . . . . . . . . . . . . . . . . . 11 D. Matrix Allocation, Deallocation and Copy . . . . . . 13 E. Display Functions . . . . . . . . . . . . . . . . . . 14 F. Matrix IO . . . . . . . . . . . . . . . . . . . . . . 15 G. Binary Matrix Functions . . . . . . . . . . . . . . . 16 H. Unary Matrix Functions . . . . . . . . . . . . . . . 17 I. Patterned Matrices . . . . . . . . . . . . . . . . . 19 J. Mathematical Functions . . . . . . . . . . . . . . . 20 V. COMPILATION AND LIMITATIONS . . . . . . . . . . . . . . . 24 VI. REVISION HISTORY . . . . . . . . . . . . . . . . . . . . 25 A. Version 1.0 . . . . . . . . . . . . . . . . . . . . 25 B. Version 1.1 . . . . . . . . . . . . . . . . . . . . 25 C. Version 1.2 . . . . . . . . . . . . . . . . . . . . 25 Mat-P 4 I. INTRODUCTION PMAT is Turbo Pascal (TP) source code for recursive matrix algebra. The program is based on another program I wrote in C to do the same thing. Both use a stack of matrices to keep track of intermediate heap allocations. Once I figured it out in C it was pretty easy to step back and do it in Pascal. The object oriented extensions in TP also helped smooth the process. This allowed a convenient scheme to allow matrices to be larger than 64K. Of course you can't overload operators in TP, so the recursion is messy. However, you can keep track of intermediate allocations on the heap using PMAT. new( a, makematrix( 1, 1 ) ); x = matequals(x, add( a,add(b,c))); This code segment has to keep track of matrix allocations on the heap, and then delete the temporary matrices. In this example, the sum of B and C is a temporary matrix which would be lost without some sort of global memory allocation tracking such as a stack. Its memory should be deleted after the equals function is called. The stack helps avoid the assignment of temporaries in recursive calls. On a more personal note, I wrote this program for fun. I started on this problem in TP 3.0 when that was the best compiler on the market (1984). I never really got a satisfactory answer, so I set it down. Then I picked it up in Turbo C, and didn't get a satisfactory answer. Then I picked it up in C++, and got a good answer in about 1990 (see YAMP on the BPROGB of Compuserve). Then I wrote it in C with good results, then TP 6.0 just for old times sake. I had forgotten how fast TP compiles. It's a real code blaster. Object oriented programming in TP 6 also helped. I guess I learned some computer science along the way too. II. FINANCIAL CONDITIONS OF USE If you like this program, send me $5.00 US (or buy a deserving beggar lunch. In either case, find a way to repay my charity.) You may distribute the package unchanged. I only plan to distribute it electronically. I retain the copyright as proof of authorship. Mat-P 5 III. BASIC ALGORITHMS A. Structures and Allocation PMAT has two important objects: vmatrix and mstack. The vmatrix object is declared in the interface vmatrixptr = ^vmatrix; vmatrix = Object r, c : integer; Function m( i, j: integer ): double; Function mm( i, j: integer ) : fp; constructor MakeMatrix( vr, vc: integer ); Destructor KillVmatrix; Procedure Garbage; Procedure Show( strng: String ); Procedure infomatrix( strng: String ); private v,vcheck: ^app; nelems : longint; onstack : boolean; Procedure purgevectors; Procedure allocvectors( rr, cc: integer ); End; The vmatrix contains integers for the dimensions of the matrices, an array of pointers to vectors of doubles, the number of elements, and a check address. Matrix allocation is based on Numerical Recipes in C . An array of addresses of vectors of doubles is allocated first. Then, an array is allocated and its first address is stored in the address vector. This method allows matrices larger than 64K and scatters the rows of the matrix throughout the heap. A row is restricted to 8190 doubles to keep it under 64K. Object oriented programming helped hide the ugly notation required for storing and retrieving elements from a vmatrix. You use the member function mm(i,j) to store elements in the matrix, and m(i,j) to retrieve elements. An example is x^.mm(i,j)^ := y^.m(i,j); where x and y are vmatrix pointers. Note the extra ^ in x^.mm(i,j)^. This returns the address in the heap for storing an element. Using x^.mm(i,j) won't work. Think of element assignment as saying "put y^.m(i,j) in the address pointed to by x^.mm(i,j)". Using m(i,j) returns a double via Mat-P 6 m := v^[i]^[j]; and using mm(i,j) returns an address via mm := @v^[i]^[j]; Using member functions protects you from having to type these unintuitive indexing operations. I also allows easy access to matrices larger than 64K. They also do range checking on the indexes. The check pointer contains the value of v upon allocating the matrix. Then it is set to zero upon deallocation of v. The check value is used in the member function Garbage(). If v does not equal vcheck, then the matrix has been deallocated, or memory may have been corrupted. The program halts if you try to use a garbage matrix (better safe than sorry!). The boolean variable, onstack, indicates that the matrix is on the stack (onstack = true) or not (onstack = false). This controls how matrices are copied. Its initial value is false. It is set to true by push(), and false by pop(). See CopySD(). PMAT expects you to use pointers to matrices, so allocation looks like this procedure junk; var x,y,z : vmatrixptr; begin new( x, makematrix( 1, 1 ) ); new( y, makematrix( 1, 1 ) ); new( z, makematrix( 1, 1 ) ); . . . { put code that uses x,y,z } . . . dispose( x, killvmatrix ); dispose( y, killvmatrix ); dispose( z, killvmatrix ); end; Call the constructor for each matrix, then make sure to dispose of the storage at the end of the procedure. Your program will develop memory leaks if you do not dispose of the vmatrix pointers. PMAT is not consistent with the MARK-RELEASE method of allocating memory. Mat-P 7 The mstack type is mstackptr = ^mstack; mstack = Object constructor InitMatrixStack; Destructor KillMatrixStack; Procedure inclevel; Procedure declevel; Function ReturnMat: vmatrixptr; Function DecReturn: vmatrixptr; Procedure push( Var a: vmatrixptr ); Procedure nrerror( strng: String ); Procedure DumpStack; private nummats, level : integer; next : mstackptr; mv : vmatrixptr; Procedure cleanstack( a: vmatrixptr ); Function pop: vmatrixptr; End; The mstack structure is a typical single linked list sort of thing. The nummats field contains the number of matrices on the stack. The level field contains the subroutine nesting level which keeps track of which temporary matrices to be disposed. As you enter and leave functions that use matrix assignment, you increment and decrement the function nesting level. The mstackptr, next, is a pointer to the next mstack object. The vmatrixptr mv is the address of a vmatrix object referenced by the stack object. B. The global stack: dispatch Mstack objects are stored in the global stack linked list with a base pointer called dispatch. Dispatch is automatically initialized in the use statement. Dereference items in dispatch using the ^. operator, such as dispatch^.level. The stack also has a tail, so if you try to deallocate it, the program stops. It has a value of nil in the field called next. C. Recursion The recursive matrix functions push mstack pointers onto the stack. The address of the matrix reference is stored in mv. You create stack elements when you call the push function. mstack elements are deallocated automatically by other functions, so you Mat-P 8 don't have to do it. Popping occurs during calls to matequals() for matrix assignment. Temporary mstack pointers are freed during matequals by the garbage collector function, CleanStack(). You do not have to pop explicitly. The only part of mstack that you have to pay attention to is the function nesting level. Garbage collection occurs during assignment. Any function that uses matrix assignment pushes matrices onto dispatch. Stack elements are deleted during garbage collection if their value of level is greater than or equal to the value of level held in dispatch. This way, each function has its own stack, but all stack elements are stored in the same linked list. You control level with the member functions, Inclevel(), Declevel(), and DecReturn(). Inclevel() and Declevel() increment or decrement level. Declevel() checks that level has not become negative. The program stops if it has. DecReturn() is used for returning matrices from functions. It decrements the nesting level, and returns the address of the matrix, v, in dispatch^.next. An example of returning a matrix pointer recursively is function f1: vmatrixptr; var x: vmatrixptr; begin dispatch^.inclevel; new( x, makematrix(1,1)); x := matequals(x, ident(5) ); dispatch^.push(x); f1 := dispatch^.decreturn; end; The stack works by manipulating pointers to vmatrix objects, rather than copies. This is faster and uses less memory. The matrix pointer is pushed onto the stack. Then the function level is decremented, and the pointer x is returned as dispatch^.next^.mv. There is no need to dispose x, since it is just a pointer stored in a global stack. The stack functions will dispose it when it is appropriate. D. Matrix Assignment Matrix assignment is accomplished using the function matequals, having prototype Function matequals(Var lop: vmatrixptr; rop: vmatrixptr) : vmatrixptr; The statement from the example above uses assignment Mat-P 9 x := matequals( x, ident(5) ); The matrix being assigned must be used in two places. The first place is in the function call where the contents of the matrix objects are manipulated. The second is on the left of the equals sign where the current function is told where the new address of x is to be found. Matequals first checks that the two inputs are not garbage. Then it checks if the two inputs point to the same 2 dimensional array. If they do, then the matrix a gets a new 2 dimensional array. The contents of rop are copied into lop. Then the intermediate stack matrices are deleted. E. Parameter Passing Parameter passing of matrices is the same as in ordinary Pascal. If you want to change the contents of a matrix that was passed as a parameter, then declare it a var parameter in the procedure. Make sure that the matrix is correctly initialized before the call. Then use matrix assignment to make the changes in the matrix. Here is an example. procedure f1( var a: vmatrixptr ); begin dispatch^.inclevel; a := matequals(a,ident(3)); dispatch^.declevel; end; procedure f2( void ); var x:vmatrixptr; begin dispatch^.inclevel; new(x,makematrix(1,1)); x^.show('before pass'); f1(x); x^.show('after pass'); dispose(x,killvmatrix); dispatch^.declevel; end; The call to f1 in f2() passes x as a var vmatrix pointer. The function call changes the contents of x from a 1x1 matrix to the 3x3 identity matrix. This example clarifies why matequals() needs to use the formal parameter name on the left of the equals sign. You have to tell the program where you put the new copy of the vmatrix. Mat-P 10 IV. FUNCTIONS There are several kinds of functions in PMAT. They fall in several categories: 0. error handling 1. array allocation or deallocation 2. stack control 3. matrix allocation, deallocation and copy 4. display 5. IO 6. binary matrix functions 7. unary matrix function 8. patterned matrices 9. mathematical functions The interface should provide a quick reference to the names and spellings. The following discussion explains the header. A. Error Handling Procedure nrerror( strng: String ); This is the error handler from Numerical Recipes in C. It is called whenever a fatal error occurs. It prints the error text and exits with a value 1. It is a member function of mstack objects. Call it as dispatch^.nrerror('something''s happening here'); Procedure Garbage; This function checks if the matrix a has been deleted or corrupted. It calls nrerror if it has. This is a vmatrix member function. Call it as x^.Garbage; Mat-P 11 B. Array Allocation or Deallocation The declarations form the structure for matrices. ap = Array[1..1] Of double; apptr = ^ap; app = Array[1..1] Of apptr; The type ap is an array of doubles. Note that the range is from 1 to 1, and the array is subscriptable and variable length. Range checking is turned off in the matrix allocation and indexing functions so that these arrays act like they do in C. Each row of the matrix is allocated separately by new(apptr). Each new apptr is stored in an subscripatble, variable length array of pointers to arrays of doubles. Procedure allocvectors( rr, cc: integer ) Arrays are allocated using the Numerical Recipes in C method for dmatrix arrays. Matrices in PMat are indexed from 1 to r, and 1 to c. procedure purgevectors. Arrays are freed using the Numerical Recipes in C method for dmatrix arrays. Neither of these functions need to be called directly since they are called by the vmatrix construction routines. They are private to the vmatrix class, so they cannot be used outside of the pmat.pas unit. C. Stack Control dispatch: mstackptr; (* stack top: global *) This declaration in the interface makes dispatch available to all units that use pmat. It is allocated automatically from the use statement. constructor InitMatrixStack; (* set up stack top and tail *) This function allocates dispatch and the stack tail. It is called by the pmat initialization routine. Procedure push( Var a: vmatrixptr ); { push a matrix onto stack } pushes a matrix onto the stack. It should be called before a f1 Mat-P 12 := dispatch^.ReturnMat; or f1 := dispatch^.DecReturn; call. The new matrix is stored in dispatch^.next^.mv. push() also sets a^.onstack to true. This is a member function of mstack. Function pop: vmatrixptr; (* free dispatch->next, return v*) pop is a private function of the matrix stack. You should not call it. It is called repetitively by CleanStack to delete temporary matrices. It deletes the matrix structure, and returns the matrix pointer mv. pop() also sets mv^.onstack to false. vmatrix *ReturnMat( void ); /* return stack top matrix */ ReturnMat() returns the value dispatch^.next^.v in functions that have not incremented the function nesting level, i.e. required matrix assignment. It does not modify the subroutine nesting level. All of the binary matrix functions push a vmatrix pointer onto the stack, and then return by the statement "f1 := dispatch^.ReturnMat;". It is a public member function of mstack. Function DecReturn: vmatrixptr; (* return stack to matrix,decrement level*) DecReturn returns the value dispatch^.next^.v in functions that have incremented the function nesting level, i.e. required matrix assignment. It decrements the subroutine nesting level. Functions that use matrix assignment, such as inv(), push matrix pointers onto the stack, and then return by the statement "f1 := dispatch^.DecReturn;". This is a public member function of mstack. Procedure inclevel; (* increment function nesting level *) Inclevel() increments the function nesting level. See the discussion of the stack and recursion. This a public member function of mstack, and is called as dispatch^.inclevel. Procedure declevel; (* decrement function nesting level *) Declevel() decrements the function nesting level. It also checks that the level has become negative. If so, it calls nrerror(). This indicates that level has been decremented more often than it has been incremented. Check the Inclevel()-Declevel() pairs, or the Inclevel()-DecReturn() pairs. See the discussion of the stack and recursion. This a public member function of mstack, and is called as dispatch^.declevel. Use this function when you use Mat-P 13 matrix assignment in a procedure, but do not want to return a matrix. Procedure cleanstack( a: vmatrixptr ); (* pop stack elements if level >= dispatch^.level *) CleanStack() is a private function of the stack. It should not be called by the user. The argument is a matrix that is not to be freed. Thus the call "CleanStack( x )" indicates that all matrices with level at or above dispatch^.level should be freed, except for matrices that have the same address as x. It calls pop while (m^.level >= dispatch^.level), and (dispatch^.next^.next != NIL). It frees all matrices not equal to x. Function matequals( Var lop: vmatrixptr; rop: vmatrixptr ) : vmatrixptr; (* copy rop into lop *) The matrix being assigned must be used in two places. The first place is in the function call where the contents of the matrix structure are manipulated. The second is on the left of the equals sign where the current function is told where the new address of x is to be found. matequals() first checks that the two inputs are not garbage. The contents of rop are copied into lop by CopySD(). Then the intermediate stack matrices are deleted by CleanStack(). D. Matrix Allocation, Deallocation and Copy constructor MakeMatrix( vr, vc: integer ); (* allocate a matrix *) MakeMatrix() constructs a matrix with vr rows and vc columns. This should be called to properly construct a matrix before using it in the matrix functions. This is a public member function of vmatrix, and should be called as new(x, makematrix(1,1)); Destructor KillVmatrix; (* deallocate a matrix *) This frees the matrix storage in m^.v, and then dispose m. It calls nrerror if m is corrupted. It is a public member function of vmatrix and should be called as dispose(m,killvmatrix); Procedure CopySD( Var source, dest: vmatrixptr ); (* copy a matrix from source to dest *) This function should be considered private to the stack function, Mat-P 14 matequals(). It is also hidden to the user in the implementation of pmat. CopySD() replaces the contents in the destination matrix by a copy of the contents in the source matrix. If source and dest are the same, then CopySD() returns immediately since you are trying to copy a matrix onto itself. CopySD() recognizes when matrices are being copied from the stack. Matrices being copied from the stack are about to be deleted by Cleanstack(), so their matrix pointer is popped into the destination matrix pointer. The destination matrix storage is released first. The indices are reset. Next, the stack matrix pointer is popped into the destination matrix pointer. Then the function returns. If the matrices point to the same dmatrix, and the source matrix is not on the stack, then a new dmatrix array is allocated for dest. If they have a different number of rows or columns, the destination 2D array is reallocated to have the same number of rows and columns. Then the dmatrix array contents of the source are copied into the destination. E. Display Functions Procedure DumpStack; This dumps the stack. It is the only function that may be called between push and the function return. It is for debugging purposes. It is a public member function of mstack and should be called as dispatch^.dumpstack; Procedure infomatrix( strng: String ); (* show dimensions *) infomatrix() displays information about the matrix. The strng is a string describing the matrix. The string, strng, is displayed, then the number of rows, columns, m, and mcheck. Use infomatrix() when the matrix is large, or when you are having problems with corrupted or deleted matrices. This is a public function of vmatrix and should be called as x^.Infomat('something is happening here'); Mat-P 15 Procedure Show( strng: String );(* display matrix *) show() displays a matrix. The string, strng is also displayed before the matrix elements. You can control the width and decimal places by calling setwid() and setdec() before calling show(). This is public function of vmatrix and a typical call is "x^.show("x");" Procedure setwid( wid: integer ); (* set display width *) Procedure setdec( decimals: integer ); (* set display decimals *) Function getwid: integer; (* get display width *) Function getdec: integer; (* get display decimals *) setwid() and setdec() store hidden global integers for the display width and decimals. show() and Writea() call getwid() and getdec() to use these variables. F. Matrix IO External matrices are stored in an ASCII character format. The first row is a title string that is at most 80 characters including the final carriage return, '\n'. The second row contains two integers separated by white space. The first integer is the number of rows, r, and the second is the number of columns, c. The next rows are the rows of the matrix. Matrix may span several text rows, but there must be r*c elements. Function reada( fid: String ): vmatrixptr; (* read ascii matrix *) Reada() reads a matrix stored in the ASCII file pointed to by fid. Reada() terminates if the title string is longer than 80 total characters, including the terminal carriage return, '\n'. Reada() should be called in a matrix assignment since it returns a matrix off of the stack. As an example, x := matequals( x, Reada('xmatrix.dat'); Procedure writea( fid: String; a: vmatrixptr; titlestr: String ); (* write ascii matrix *) Writea() stores a vmatrix, a, in file, fid, using the PMat ASCII file format. Data is written to fid using the current print display values found in setdec() and setwid(). The title string Mat-P 16 is truncated to 80 characters (including '\n') if it is longer than 80 characters. G. Binary Matrix Functions The binary matrix functions are recursive in that they push matrices onto the stack, and do not manipulate the function nesting level. They all push a vmatrix pointer onto the stack, and return a vmatrix pointer using a call to ReturnMat(). The incoming matrices are checked by calling Garbage() before doing any calculations. A typical use of these functions is (* d= a+(b-c) ; *) d = equals( d, add(a,sub(b,c))); The call to equals deletes the temporary difference of c from b. Scalars may also be used in binary operations. Each of the binary matrix functions have "sc" appended to the left(right) if a scalar is to be added to a matrix on the left(right). function add( a, b: vmatrixptr ): vmatrixptr; (* addition a+b *) function scadd( a: double; b: vmatrixptr): vmatrixptr; function addsc( b: vmatrixptr; a: double); vmatrixptr; Matrix addition, conditions : a^.r == b^.r and a^.c == b^.c The scalar additions add a double to a matrix. Function sub( A, B: vmatrixptr ):vmatrixptr; (* subtraction a-b*) Function scsub( a: double; b: vmatrixptr): vmatrixptr; Function subsc( b: vmatrixptr; a: double); vmatrixptr; Matrix subtraction, conditions : a^.r = b^.r and a^.c = b^.c The scalar subtractions perform the subraction in the indicated order: scsub(i,j) = A - b(i,j), subsc(i,j)= b(i,j)-A. Function Mult( A, B: vmatrixptr ):vmatrixptr;(*matrix mult a*b*) Caley multiplication, inner products are formed using extended doubles, then truncated back to doubles when stored in a matrix. conditions : a^.c = b^.r. Mat-P 17 Function emult( a, b: vmatrixptr ):vmatrixptr; (* elementwise mult a#b *) function scmult( a: double; b: vmatrixptr): vmatrixptr; function multsc( b: vmatrixptr; a: double); vmatrixptr; Elementwise multiplication, conditions : a^.r = b^.r and a^.c = b^.c The scalar multiplications multiply the elements of b by a. Function divis(a,b:vmatrixptr): vmatrixptr; function scdivis( a: double; b: vmatrixptr): vmatrixptr; function divissc( b: vmatrixptr; a: double); vmatrixptr; Elementwise division, conditions : a^.r = b^.r and a^.c = b^.c and for all i,j b^.m(i,j) <> 0 The scalar divisions perform division in the indicated order: scdivis(i,j) = A/b(i,j), divissc(i,j)= b(i,j)/A. H. Unary Matrix Functions The unary matrix functions, neg(), tran(), and inv(), are recursive in that they push matrices onto the stack, and do not manipulate the function nesting level. neg(), and tran() push a vmatrix pointer onto the stack, and return a vmatrix pointer using a call to ReturnMat(). inv() uses matrix assignment so it increments the nesting level, pushes the return matrix, and returns by a call to DecReturn(). The incoming matrices are checked by calling Garbage() before doing any calculations. A typical use of these functions is Mat-P 18 function regression: vmatrixptr; var x,y,xpx,xpy,beta : vmatrixptr; begin (* beta = inv(x'x)x'y *) dispatch^.Inclevel; new(x, makematrix(1,1)); new(y, makematrix(1,1)); new(xpx, makematrix(1,1)); new(xpy, makematrix(1,1)); new(beta, makematrix(1,1)); x := matequals( x, Reada( "x.dat") ); y := matequals( y, Reada( "y.dat") ); xpx := matequals( xpx, mult( tran(x), x )); xpy := matequals( xpy, mult( tran(x), y )); beta := matequals( beta, mult( inv(xpx), xpy)); dispose( x, killvmatrix); dispose( y, killvmatrix); dispose( xpx, killvmatrix); dispose( xpy, killvmatrix); dispatch^.push( beta ); regression := dispatch^.DecReturn; end; function sweepreg: vmatrixptr; var xy, xypxy, beta : vmatrixptr; begin (* beta using sweep *) dispatch^.Inclevel; new(xy, makematrix(1,1)); new(xypxy, makematrix(1,1)); new(beta, makematrix(1,1)); xy := matequals( xy, Reada( "xy.dat") ); xypxy := matequals( xypxy, mult( tran(xy), xy )); xypxy := equals( xypxy, sweep( xypxy, 1, xy^.c - 1) ); beta := matequals( beta, submat( xypxy, 1, xy^.c -1, xy^.c, xy^.c) ); dispose( xy, killvmatrix); dispose( xypxy, killvmatrix); dispatch^.push( beta ); sweepreg := dispatch^.DecReturn; end; Function neg( A: vmatrixptr ): vmatrixptr;(* negation *) Changes the sign of all elements. Function tran( a: vmatrixptr ): vmatrixptr;(* transpose *) Matrix transpose: c^.m(i,j) = a^.m(j,i) Mat-P 19 Function sweep( A: vmatrixptr; k1, k2: integer ) : vmatrixptr; Sweep cols k1...k2 along the diagonal and in place in a. See the example for how to use sweep for regression. Input matrix must be square. Row and columns having pivot elements less than 10E-8 are set to zero. This indicates a colinearity with at least one column that has already been swept out. Function inv( Amat: vmatrixptr ): vmatrixptr; (* inversion by sweep *) Calls sweep for inversion: x = matequals( x, sweep(x, 1, x^.c )); This is a Gaussian elimination inverter without column pivoting. It fails to accurately invert Hilbert matrices of dimension 8 or more, which is reasonable for inversion by Gaussian elimination. I. Patterned Matrices The patterned matrix functions are recursive. They push matrices onto the stack and return by calling ReturnMat. They should be used in conjunction with matrix assignment. Function submat( a: vmatrixptr; lr, r, lc, c: integer ): vmatrixptr; (*submatrix*) Extract a submatrix of a, using rows lr to r, and columns lc to c. All integer arguments must be in the range of the dimension of a. Function ident( n: integer ): vmatrixptr; (*identity matrix*) Construct an n dimensional Identity matrix Function fill( rr, cc: integer; d: double ):vmatrixptr; (* matrix of constants d *) Construct a rr x cc matrix of constants, d. Function cv( a, b: vmatrixptr ): vmatrixptr; (* vertical concatenation *) Vertically concatenate b to a. a is on top, b is on bottom, and they must have the same number of columns. Mat-P 20 Function ch( a, b: vmatrixptr ): vmatrixptr; (* horizontal concatenation *) Horizontally concatenate b to a. a is on the left, b is on the right, and they must have the same number of rows. vmatrix *vecdiag( vmatrix *a ); (* place vector on diagonal *) Place elements of a rx1 or 1xc matrix on the diagonal of a square matrix. J. Mathematical Functions The mathematical function module for YAMP was translated for PMAT. These are some of the many linear algebra functions that are used. function MSort(a :vmatrixptr; col, order: integer): vmatrixptr; Sort the rows of a matrix by a column if order is not equal to zero. If order is zero, then sort the columns of a matrix by a row. function Mexp(ROp: vmatrixptr): vmatrixptr; function Mabs(ROp: vmatrixptr): vmatrixptr; function Mlog(ROp :vmatrixptr): vmatrixptr; function Msqrt(ROp: vmatrixptr): vmatrixptr; These functions operate on the elements of a matrix. They halt if an argument is out of range. The take exp(a^.m(i,j)), abs(a^.m(i,j)), ln( a^.m(i,j)), and sqrt(a^.m(i,j)). You might consider adding the trig functions too. function Trace(ROp: vmatrixptr): double; Takes the sum of the diagonal elements of a square matrix. function Helm (n: integer): vmatrixptr; Returns the Helmert matrix of dimension n. function Index( lower, upper, step : integer): vmatrixptr; Mat-P 21 Returns an index matrix with elements between lower and upper, and in increments of step. If lower is less than upper, then step must be positive. If lower is greater than upper, then step must be negative. function Kron(a,b:vmatrixptr): vmatrixptr; Returns the Kroniker product of a and b. The blocks of the Kroniker product are a^.m(i,j)*b. Mat-P 22 function Det(Datain : vmatrixptr): double; Returns the determinant of Datain, which must be a square matrix. function Cholesky(ROp: vmatrixptr): vmatrixptr; Returns the Cholesky decomposition of a square symmetric matrix: ROp = S'S, where S is an upper triangular matrix. Cholesky stops if ROp is singular. This routine does not use pivoting. procedure QR(var ROp, Q, R: vmatrixptr; makeq: boolean); This procedure produces the QR decomposition of an rxc matrix: ROp=QR, where Q is an rxc matrix of orthogonal columns, and R is an upper triangular cxc matrix. This routine stops if a zero is detected along the diagonal of R. The variable makeq determines if Q is produced. The procedure calculates Q = ROp*Inv(R). If the calculation of Q is not required, or inversion of R is prohibative, then set makeq = FALSE to avoid the inversion of R. function Svd(var A, U, S, V : vmatrixptr; makeu, makev: boolean) : integer; This function returns an integer indicating that the singular value decomposition has failed. The singular value decomposition is calculated as A = UDiag(S)V' where U and V are orthogonal matrices and S is a column vector containing the singular values of A. The singular values of A are the eigenvalues of A'A. If you wish to calculate U or V, then set makeu = TRUE or makev = TRUE respectively. function Ginv( ROp : vmatrixptr): vmatrixptr; This function returns the generalized inverse from the singular value decomposition of ROp=UDiag(S)V'. The generalized inverse of ROp is Ginv(ROp) = VDiag(S+)U', where S+ is a diagonal matrix of the reciprocals of S. If S has a element of zero, the corresponding element of S+ is also zero. function Fft(ROp : vmatrixptr; INZEE : boolean): vmatrixptr; This function returns the fast Fourier transform of an input matrix ROp. The length of ROp may be any value. ROp may have 1 or 2 columns. If it has 1 column, the vector is assumed to be a real vector. If the matrix has 2 columns, the matrix is assumed to consist of complex numbers. The first column contains the real part, and the second column contains the imaginary part. The Mat-P 23 forward transform is calculated if INZEE is TRUE, and the inverse transform is calculated if INZEE is FALSE. function Vec( ROp : vmatrixptr): vmatrixptr; This functions stacks the columns of a matrix into a single column vector. Thus if ROp is rxc, then Vec(ROp) has r*c rows and 1 column. The c blocks of Vec(ROp) are columns of ROp. function Diag(ROp : vmatrixptr): vmatrixptr; This function has two purposes. If the input matrix is a column or row vector, then Diag places the elements along the diagonal elements of square matrix of zeros. If the input matrix is square, then all of the off-diagonal elements are set to zero. function Shape(ROp: vmatrixptr; nrows: integer): vmatrixptr; This functions reshapes a matrix to have nrows rows. If ROp is has r*c elements, then nrows must divide r*c without a remainder. type margins = (ALL,ROWS,COLUMNS); This enumerated type controls the actions of the next 5 functions. If you use ALL, then the functions operate over all elements of the input matrix. If you use ROWS, then the functions operate over the rows of the input. If you use COLUMNS, then the functions operate over the columns of the matrix. function Sum( ROp : vmatrixptr; marg : margins ) : vmatrixptr; This functions returns sum of the elements in a matrix. If marg=ALL, then the returned matrix is 1x1, containing the sum of all elements. If marg=ROWS then the return matrix has the same number of columns as ROp, and 1 row. Each column contains the sum of the row elements. If marg=COLUMNS then the return matrix has the same number of rows as ROp, and 1 column. Each row cantains the sum of the column elements. function Sumsq( ROp : vmatrixptr; marg: margins): vmatrixptr; This function is similar to Sum, except that it returns the sum of squared elements. Mat-P 24 function Cusum( ROp : vmatrixptr): vmatrixptr; This function returns the cumulative sum across rows of matrix. function Mmax( ROp: vmatrixptr; marg : margins): vmatrixptr; function Mmin( ROp: vmatrixptr; marg : margins): vmatrixptr; These functions find the max and min elements of a matrix, and the position of the elements. If marg=ALL, then the functions return a 3x1 matrix. The first element is the row index of the max, the second element is the column index of the element, and the third element is the max(or min). If marg=ROWS, then the returned matrix is a 3xc matrix where the first row is the row index of the max(min) element in a column, the second row is the column index of the max(min) element, and the third row are the max(min) for the column. Similarly, if marg=COLUMNS, then a rx3 matrix is returned. Procedure Crow(var ROp: vmatrixptr; row : integer; c : double); Procedure Srow(var ROp : vmatrixptr; row1, row2 : integer); Procedure Lrow(var ROp : vmatrixptr; row1, row2 : integer; c : double); Procedure Ccol(var ROp : vmatrixptr; col : integer; c : double); Procedure Scol(var ROp : vmatrixptr; col1, col2 : integer); Procedure Lcol(var ROp : vmatrixptr; col1, col2 : integer; c : double); These functions perform elementary row and column operations on matrices. They change the input matrices. The first letter indicates the operation. The next 3 letters indicate if the operation is for rows or columns. For Crow, the row indexed by row is multiplied by c. For Srow, row1 and row2 are swapped. For Lrow, row1 is multiplied by c, and added to row2. Ccol, Scol, and Lcol perform similar operations on columns. V. COMPILATION AND LIMITATIONS Programs using PMAT must include pmat in the uses statement. Also put pmatop in the uses statement if you use code from PMATOP.PAS. Compile and link the program. See mattest2.pas or testop.pas for example files using PMat or PMatop. Also see makefile for a Borland make file for using the command line compiler. Mat-P 25 PMat does have some limitations. It is limited by the available memory in the machine, so large matrices will exhaust the heap, and the program will stop. This is bothersome on a PC without a DOS extender. The BP7 protected mode for DOS programs allows up to 16 megs of memory, and matrices may be larger than 640K. Note that a 1000x1000 of double takes 8 megs, so PMat may not be suitable for really large matrices. Another limitation is its readability. Using function calls rather than overloaded algebraic operators muddles things. There is no other obvious way to do this in TP though. One trick I use to help me with this is to count open-close parenthesis pairs. Some programming editors help this too by blink the other member of the pair. Another trick is to keep recursive expressions on a single line. This helps with storage also, since the stack storage is cleaned up after each call to equals. If you have memory overflow problems, try splitting the recursive expressions into several parts by paying attention to where the big matrices are created. VI. REVISION HISTORY A. Version 1.0 This is the original version. It is a copy of my matrix program C-Mat. I may extend PMAT to have most of the capability of YAMP, my C++ program. This would involve a virtual memory manager, a graphics module, and larger set of mathematical functions. I don't know if it is worth the effort though. B. Version 1.1 This version adds BP 7.0 compatability for protected mode DOS applications. The version 7 DPMI extender eliminates the need of a virtual memory manager. You need a 286 or 386 to use this feature. C. Version 1.2 Added a new module for the math functions from YAMP. This includes matrix sorts, elementwise functions, trace, helmert matrix, index matrix, Kroniker products, determinants, Cholesky decomposition, QR decomposition, singular valued decomposition, generalized inverse, fast Fourier transform, vec, diagonal, shape, sums, sums of squares, cumulative sums, max, min, and elementary row-column operations. Added arithmetic operations with scalars: scadd, addsc, etc.